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Abstract 

Strictly speaking the laws of the conventional 
Statistical Physics, based on the Equipartition 
Postulate [T] and Ergodicity Hypothesis [2) . apply 
only in the presence of a heat bath. Until recently 
this restriction was not important for real phys¬ 
ical systems: a weak coupling with the bath was 
believed to be sufficient. However, the progress 
in both quantum gases and solid state coherent 
quantum devices demonstrates that the coupling 
to the bath can be reduced dramatically. To de¬ 
scribe such systems properly one should revisit 
the very foundations of the Statistical Mechanics. 
We examine this general problem for the case of 
the Josephson junction chain and show that it dis¬ 
plays a novel high temperature non-ergodic phase 
with finite resistance. With further increase of the 
temperature the system undergoes a transition to 
the fully localized state characterized by infinite 
resistance and exponentially long relaxation. 

The remarkable feature of the closed quantum sys¬ 
tems is the appearance of Many-Body Localization 
(MBL)[3]: under certain conditions the states of a 
many-body system are localized in the Hilbert space 
resembling the celebrated Anderson Localization [3] 
of single particle states in a random potential. MBL 
implies that macroscopic states of an isolated system 
depend on the initial conditions i.e. the time averag¬ 
ing does not result in equipartition distribution and 


the entropy never reaches its thermodynamic value. 
Variation of macroscopic parameters, e.g. tempera¬ 
ture, can delocalize the many body state. However, 
the delocalization does not imply the recovery of the 
equipartition. Such a non-ergodic behavior in iso¬ 
lated physical systems is the subject of this Letter. 

We argue that regular Josephson junction arrays 
(JJA) under the conditions that are feasible to imple¬ 
ment and control experimentally demonstrate both 
MBL and non-ergodic behavior. A great advantage 
of the Josephson circuits is the possibility to disentan¬ 
gle them from the environment as was demonstrated 
by the quantum information devices. [S] At low tem¬ 
peratures the conductivity a of JJA is finite (below 
we call such behavior metallic), while as T —>■ 0 JJA 
becomes either a superconductor (cr —)■ oo) or an insu¬ 
lator (tr —>• 0)[S1[7]. We predict that at a critical tem¬ 
perature T = Tc JJA undergoes a true phase transi¬ 
tion into a MBL insulator (ct = 0 for T > Tc). Re¬ 
markably already in the metallic state JJA becomes 
nonergodic and can not be properly described by the 
conventional Statistical Mechanics. 

JJA is characterized by the set of phases {cfn} and 
charges{( 7 i} of the superconducting islands, ipi and qi 
for each i are canonically conjugated. The Hamilto¬ 
nian H is the combination of the charging energies 
of the islands with the Josephson coupling energies. 
Assuming that the ground capacitance of the islands 
dominates their mutual capacitances (this assump- 
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tion is not crucial for the qualitative conclusions) we 
can write H as 




+Ej{l- cos( 


( 1 ) 


The ground state of the model 0 is determined 
by the ratio of the Josephson and charging energies, 
Ej /Ec that controls the strength of quantum fluctu¬ 
ations: JJA is an insulator at Ej/Eq < ij and a su¬ 
perconductor at Ej/Ec > ?7 [ni[7] with rj « 0.63 (see 
Supplemental materials). The quantum transition 
at Ej/Ec = rj belongs to the Berezinsky-Kosterlitz- 
Thouless universality class [8]. Away from the ground 
state in addition to Ej/Ec there appears dimension¬ 
less parameter U/Ej where U is the energy per su¬ 
perconducting island {U = T in the thermodynamic 
equilibrium at T <C Ej). 

The main qualitative finding of this paper is the ap¬ 
pearance of a non-ergodic and highly resistive “bad 
metal” phase at high temperatures, T/Ej > 1, which 
at T > Tc « E'^j/Ec undergoes the transition to the 
MBL insulator. In contrast to the T = 0 behavior, 
these results are robust, e.g. are insensitive to the 
presence of static random charges. The full phase 
diagram in the variables Ej/Ec , T/Ej is shown in 
Fig.[Tl We confirmed numerically that the bad metal 
persists in the classical limit although Tc —> oo; it is 
characterized by the exponential growth of the resis¬ 
tance with T and violation of thermodynamic identi¬ 
ties. We support these findings by semi-quantitative 
theoretical arguments. Finally, we present the results 
of numerical diagonalization and tDMRG (time Den¬ 
sity Matrix Renormalization Group jS]) of quantum 
systems that demonstrate both the non-ergodic bad 
metal and the MBL insulator. 

It is natural to compare the non-ergodic state of 
JJA with a conventional glass characterized by in¬ 
finitely many metastable states. The glass entropy 
does not vanish at T = 0, i. e. when heated 
from T = 0 to the melting temperature the glass re¬ 
leases less entropy than the crystal (Kauzmann para¬ 
dox [10]). Similarly to glasses JJA demonstrates 
non-ergodic behavior in both quantum and classical 
regimes. However the ergodicity violation emerges 
as high rather than low temperatures transforming 



Fig. 1 : Phase diagram of one dimensional Joseph¬ 
son junction array. The MBL phase transi¬ 
tion separates the non-ergodic bad metal with 
exponentially large but finite resistance from 
the insulator with infinite resistance. Cooling 
the non-ergodic bad metal transforms it into 
a good ergodic metal. The points show ap¬ 
proximate positions of the effective T/Ej for 
the quantum problem with a finite number of 
charging states. The red stars indicate insula¬ 
tor, blue circles bad metal, and squares good 
metal. 


Kauzmann paradox into an apparent temperature di¬ 
vergence (see below). 

Qualitative arguments for MBL transition. In a 
highly excited state U Ej the charging energy 
dominates: Ecq^ ^ U ^ Ej. Accordingly, the value 
of the charge, \qi\, and charge difference on neigh¬ 
boring sites, 5qi = \qi — qi+i\ are of the order of 
q^5q^ The energy cost of a unit charge 

transfer between two sites SE ~ \/UEc exceeds the 
matrix element of the charge transfer, Ej/2, as long 
asU ^ Umbl = E'^j/Ec- According to |11I3] the sys¬ 
tem is a non-ergodic MBL insulator under this condi¬ 
tion. Thus, we expect the transition to MBL phase at 
Tc/Ej cx Ej/Ec with the numerical prefactor close 
to unity (see supplemental materials). 

If Ej/Ec 1 and U ^ T Ej the conductivity 
limited by thermally activated phase slips is exponen¬ 
tially large, a ^ exjo(E t /T) [TTl fT^ . As we show be- 
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low, at T ^ Ej in the metallic phase conductivity is 
exponentially small, a ^ exp(—T/i?j), even far from 
the transition, Ej T Ej/Eq- The resistance 
in this state can exceed Rq = /i/(2e)^dramatically 
and still display the “metallic” temperature behavior 
{dR/dT > 0). 

Classical regime is realized at Ec —>■ 0 for fixed Ej 
and T. One can express the charge of an island q 
through the dimensionless time r = y/EjEct as q = 
^JEJ/Ecd4>/dT. Since q ^ ^Ej/Ec > 1 . one can 
neglect the charge quantization and use the equations 
of motion 


dr'^ 


= sm{(j)i+i - (j)i) + sin((;!)j_i - ^j). 


( 2 ) 


Here i = and the boundary conditions are 

00 = 4’l+i = 0. We solve these equations for the var¬ 
ious initial conditions corresponding to a given total 
energy and compute the energy Us contained in a 
part of the whole chain of the length 1 <C ? T as a 
function of r. 

The ergodicity implies familiar thermodynamic 
identities, e.g. 

[{Ul)-{Usf)/T^ = d{Us)/dT (3) 


This relation between the average energy of the sub¬ 
system, (Us) and its second moment (C/|) turns out 
to be invalid for a bad metal. To demonstrate this we 
evaluated the average energy per site in this subsys¬ 
tem, u = {Us)t/[EjI) and the temporal fluctuations 


of this energy, Wr = - {U/ {E^l). Here 

(.. .)^ and the bar denote correspondingly averaging 
over the time Q and over the ensemble of the initial 
conditions. 

From Q it follows that u = T/{2Ej) at T ^ Ej 
{u (T)-function is evaluated for arbitrary T /Ej in 
supplemental materials). One can thus rewrite @ 
as 


w = 


T 2 du 


2u^ 


( 4 ) 


Results of the numerical solution of ([^ are com¬ 
pared with ([^ in Fig. For any given evolution 


^ Given the evolution time Tav (• • •)t- averaging means aver¬ 
aging over the time interval Tev after initial evolution for time 
'Eev • 



Fig. 2: Energy fluctuations as a function of average 
energy per island for different evolution times 
{Tey = 1, 2, 4, 8,16 X 10^ ) in a small sub¬ 
system; the dashed line is the extrapolation 
to infinite times as explained in the text. A 
single point (red) obtained by direct compu¬ 
tations up to Tev = 2 10 ® at which time scales 
the time dependence practically disappears 
for L < 100. The upper (red) curve cor¬ 
responds to the thermodynamic identity (§. 
The insert shows the effective temperature de¬ 
fined by the energy fluctuations determined 
at different time scales: Tev = 2, 4, 8 x 10"*^ 
(Ti — T 3 respectively) and by their extrapo¬ 
lation to infinite times (T*), their comparison 
with the temperature expected in thermody¬ 
namic equilibrium, Tth- 


time, Tev , the computed Wr (u)- dependence saturates 
instead of increasing as Q. At a fixed u, Wr (u) 
increases with time extremely slowly. Below we ar¬ 
gue that Wr (u) does not reach Its thermodynamic 
value even at r —>■ 00 . Violation of the thermody¬ 
namic identity implies that temperature is ill defined, 
so the average energy u rather than T is the proper 
control parameter. The effective temperature defined 
as T*(tt) = Ej j Jq du/w{u), is shown in the insert to 

Fig-d 

Qualitative interpretation. Large dispersion of 
charges on adjacent islands, z, i -I- 1 at u 3> 1 im¬ 
plies quick change of the phase differences, (pi — 0 i+i, 
with time. Typical current between the two neigh¬ 
bors Ajsin( 0 i — 0 i+i) time-averages almost to zero. 
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However, accidentally the frequencies, uji = d4)i/dt, 
get close. In the classical limit the difference uji—uji+i 
can be arbitrary small. Such pair of islands is char¬ 
acterized by one periodic in time phase difference. 
Contrarily, three consecutive islands with close fre¬ 
quencies \u!i — Wi+il ^ \uji — Wi-il < 1 experience 
chaotic dynamics that contains arbitrary small fre¬ 
quencies, similarly to work [13) . For uncorrelated 
frequencies with variances u ^ 1 , a triad of islands 
(j — 1 , i, i-1-1) is chaotic with the probability 1/u ^ 1, 
i.e. such triads are separated by large quiet regions 
of a typical size rj ~ u ^ 1. The low frequency noise 
generated by a chaotic triad decreases exponentially 
deep inside a quiet region. Provided that lu uji+m 
, m = 0 ... r the superconducting order parameter 
Zi{uj) = f dTexp(i(^i + icur) at site i satisfies the re¬ 
cursion relation 

r-l ^ 

Z,+rH = Zi(uj) ^- ( 5 ) 

m—0 

which implies the log-normal distribution for the re¬ 
sistances Rj,j+r of quiet regions (see Supplemental 
materials): 

(In^ (i?/i?t)) = In Rt (6) 

lni?t = (Ini?) = itln(u) (7) 


where Rt is the typical resistance of a quiet region. 
The resistance of the whole array is the sum of the 
resistances of the quiet regions. The mean number of 
these regions in the chain equals to N = L/rt ^ 1, 
its fluctuations being negligible. For the log normal 
distribution the average resistance of a quiet region, 
(i?), is given by(i?) = For the resistivity, p, we 

thus have p = N (R) jL = R^^^Ivt- According to Q 

3 

In prt « - (7 -I- In u) rt « u(ln u + ^) ( 8 ) 


where 7 = 0.577 is Euler constant. 

In order to determine the current caused by voltage 
V across the chain we solved the equations (§ with 
modified boundary conditions, (/>o = 0 , — Vt. 

The results confirm the prediction (H), see Fig. § 
The range of the resistances set by realistic com¬ 
putation time is too small to detect the logarith¬ 
mic factor in (8[7), however, a relatively large slope. 



Fig. 3; Resistivity as a function of the internal en¬ 
ergy. At high energies resistivity is expo¬ 
nentially large due to large regions of almost 
frozen charges (see text). Insert: at low 
energies (temperatures) exponentially large 
conductivity is limited by exponentially rare 
phase slips. High u points require computa¬ 
tion times Tev ~ 10 ®. 


d\np/du « 3.0 at it = 3.5 is consistent with ([^ that 
gives dlnp/du « 2.5 -I- 1.5 In it — l/it. 

The qualitative picture of triads separated by log- 
normally distributed resistances of silent regions al¬ 
lows one to understand the long time relaxation of 
Wt{u) in the subsystem of length I (see Fig. [^. Each 
resistance can be viewed as a barrier with a tunnel¬ 
ing rate ~ 1/R. For a given time t the barriers with 
T R can be considered impenetrable, whilst the 
barriers with t ^ R can be neglected. As a result, 
the barriers with R> t break the system into essen¬ 
tially independent quasiequilibrium regions (QER) of 
the typical size 


Ir 


n 

, = exp 

V27r In Rt 


In {r/Rt) 
2 In i?t 


(9) 


li I ^ rt and r < exp A/ln(l/rt) In Rr the subsys¬ 
tem contains It/I ^ 1 QER, so that w oc. It/I- At 
longer times, It ^ I, the subsystem is in equilibrium 
with a particular QER and w (x I/It- The full depen¬ 
dence on time can be interpolated as 


Wo 


Wt = 


1 -|-/3exp (—Qfln^ (t/to)) 


( 10 ) 
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Fig. 4: Time dependence of energy fluctuations 
for the subsystems of different sizes I = 
10, 20, 40, 80 (left to right) and their fits to 
log-normal law ( [l0| . The length of the full 
system was L = 51. The best fit shown here 
corresponds to the value a = 0.05, (3/1 = 
1.0 — 1.5 and Inrg « 3.0 — 4.5. The values for 
larger sizes {jS/l « 1.0 and Inro « 4.5 agree 
very well with the ones expected for Rt ob¬ 
tained in the computation of the resistance at 
u = 3.5 dni?( « 5.0. This yields Tt « 5.0 , 
P/l « 1.1 and Inro ~ 5.0. 


where a = (21ni?t) /3 = \pKja{J,lrt) and tq = Rt- 

The numerical simulations confirm that the energy 
variance w relaxes in agreement with (10) as shown 
in Fig. The best fit to the equation ( (l0| yields 
parameters close to the expected, lni?( = \-a.p(u)rt, 
rt ^ u. Extrapolation gives Woo(oo) « 10.0 of Woo{l) 
to large I, which is significantly smaller than thermo¬ 
dynamic value WTh = 14.0 at m = 3.5. 

Another test of the ergodicity follows from the 
fluctuation-dissipation theorem (FDT) that relates 
conductivity and current fluctuations. In the low 
frequency limit the noise power spectrum is S{uj —> 
0 ) = 2Teff/R where T^ff is the effective temperature 
(see Fig. [^, which we extracted from the numerical 
data. We found that Tg// > Tth for m > 1, where 
Tth is thermodynamic temperature. In particular, 
Teff « l.dTxh for u = 3.5, which is close to Ti.(m) 
shown on the inset to Fig. 

Quantum behavior. In contrast to a classical limit 
Ec —>■ 0 in the quantum regime Ec > 0 we ex¬ 



Fig. 5; Spectrum of current noise fluctuations, 
5'(a;) for different internal energies, U = 
2.5, 3.0, 3.55. The dashed line shows S{uj) ~ 
dependence. Insert shows the effective 
temperature determined from FDT relation. 


pect at T = Tc a MBL phase transition between 
two non-ergodic states: the insulator {p = oo) and 
a bad metal {p < oo). For Ej ^ Ec the bad metal 
can be described classically at T ^ Tc ~ E'j/Ec. 
Our previous discussion suggests that the bad metal 
is non-ergodic in a broad range of the parameters, 
T/Ej and Ej/Eq- To verify this conjecture numer¬ 
ically we reduced the Hilbert space of the model 0 
to a finite number of charging states at each site, 
Qi = 0,±1,±2 (RHS model). We analyzed the time 
evolution of entanglement entropy, 5'{'I'} of the left 
half of the system. The entropy was averaged over 
the initial states from the ensemble of product states 
in the charge basis, St{L) =< S'j'I'} that corre¬ 
spond to zero total charge. As a result, we obtained 
the Gibbs entropy at T = oo (all states have the same 
weight exp{—H/T) = I) 

The insert to Fig. shows the time dependence 
of the entropy at Ej/Ec = 0.3 which corresponds 
to the bad metal regime (see below). A slow satura¬ 
tion of the entropy follows its quick initial increase. 
It is crucial that the saturation constant, Soo (L) is 
significantly less than its maximal value, SxhiL) = 
Lin5 expected at T = oo equilibrium. Furthermore, 
dSao(,L)/dL < In5, indicating that Sao — Sth is ex¬ 
tensive and the system is essentially non-ergodic. 

Fig. [^presents Sao as a function of Ej/Eq. Note 
that Sao is measurably less than Sxh for Ej/Eq < 
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Fig. 6: Enropy per spin in thermodynamic limit ex¬ 
trapolated to t = 00 . The insert shows 
S{t) dependence in the bad metal regime at 
Ej/Ec = 0.3 that shows slow relaxation of 
the entropy to its saturation value. 

0.6 — 1. For Ej/Ec ^ 1, the entropy saturation is 
quick and the accuracy of the simulations does not 
allow us to distinguish Sao from Sth (see Supple¬ 
menting material). We thus are unable to conclude 
if the system is truly ergodic or weakly non-ergodic 
at Ej/Ec ^ 1. The former behavior would imply a 
genuine phase transition between bad and good met¬ 
als, while the latter corresponds to a crossover. 

Deep in the insulator the time dependence of the 
entropy is extremely slow, roughly linear in In t in a 
wide time interval (see Supplemental materials). This 
resembles the results of the works [H US] for the 
conventional disordered insulators. The extremely 
long relaxation times can be attributed to rare pairs 
of almost degenerate states localized within different 
halves of the system. The exponential decay with dis¬ 
tance of the tunneling amplitude that entangles them 
leads to the exponentially slow relaxation. 

In order to locate the MBL transition we analyzed 
the time dependence of the charge fluctuations. In a 
metal the charge fluctuations relaxation rate depends 
weakly (as a power law) on the sample size in contrast 
to the exponential dependence in the insulator. Com¬ 
paring the dependences of the rates on the system size 
for different Ej/Ec (see Fig. 0 ) we see that the tran¬ 


sition happens in the interval 0.05 < Ej/Ec < 0.3. 

The variances of the charge in the RHS model at 
T = oo and the problem 0 at finite T coincide 
at T = 2Ec- Thus, we expect that the results of 
the quantum simulations describe the behavior of 
the model Q at T/Ej ~ Ec/Ej yielding the hy¬ 
perbola shown in Fig. The MBL transition at 
Ej/Ec ^ 0.2 in T = oo RHS model corresponds to 
the transition temperature Tc ~ lOEj in model 0. 
The transition line shown in Fig. [^is a natural con¬ 
nection of this point with « Ej/Ec asymptotic at 
Ej/Ec ^ 1 discussed above. The maximum of the 
transition temperature is natural because the quan¬ 
tum fluctuations are largest at Ej ~ Ec- 

Possible experimental realization. MBL and the 
violation of the ergodicity can be observed only at 
sufficiently low temperatures when one can neglect 
the effects of thermally excited quasiparticles which 
form the environment to model Q . This limits tem¬ 
peratures to r < O.IA where A is superconduct¬ 
ing gap. In order to explore the phase diagram 
one has to vary both T/Ej and Ej/E q in the in¬ 
tervals 1 < T/Ej < 5 and 0.1 < Ej/Ec < 5. 
The former condition can be satisfied if each junc¬ 
tion is implemented as a SQUID loop with individ¬ 
ual Josephson energy E^p ^ T^ax = O.IA so that 
Ej = 2 cos(e$/ 7 rIi)Fj where $ is flux through the 
loop. The latter condition requires enhancing ground 
capacitance of each island which should exceed the 
capacitance of the junctions in order for the model 
Q to be relevant. Realistic measurements of such 
array include resistance R{T^ Ej) and current noise. 
Here we predict a fast growth with temperature and 
divergence at of the resistivity and violation of 
FDT. 

In conclusion, we presented strong numerical 
evidences for the MBL transition and its semi- 
quantitative description in a regular, disorder fre^ 
Josephson chain. Probably the most exciting find¬ 
ing is the intermediate non-ergodic conducting phase 
(bad metal) between the MBL insulator and good 
ergodic metal. This phase distinguishes Joseph- 

^ When this work was finished we learnt about the numerical 
studies [Tfil I17| that reports MBL in disorder free ID systems, 
in cotrast the work m reports the absence of MBL in such 
systems. All these works studied models different from ours. 
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Fig. 7: Charge relaxation in a bad metal {Ej/Ec = 
0.3). The characteristic time scale, ter, can 
be defined by the extrapolated crossing point 
of (n^) with t-axis. The crossing point ter 
shows dramatically different behavior for the 
insulator and metal as shown by insert: in 
the good metal this time is very short, in the 
bad metal it is dramatically longer but does 
not increase with size whilst in insulator it is 
extremely long and grows with size. 


son junction chain from the spinl/2 Heisenberg one¬ 
dimensional model in the random field, na uni EH, 
where the ergodicity is believed to be violated only 
in the MBL regime [51]. The ergodicity violation in 
a wide range of parameters is in contrast to the one 
particle Anderson localization in a finite-dimensional 
space where the non-ergodic behavior is limited to the 
critical point. However, the signatures of the “non- 
Gibbs regime” in the discrete nonlinear Schrodinger 
equation [55] |53| , the appearance of the non-ergodic 
states in a single particle problem on the Bethe lat¬ 
tice as well as the dominance of single path relaxation 
close to the critical point of superconductor-insulator 
transitions on Bethe lattice [531 Hi] rnake us believe 
that the intermediate non-ergodic phase is a generic 
property of MBL transition rather than an exception. 
Further work is needed to establish the domain of the 
applicability of these results as well as the nature of 
the transition between bad and good metals. 

Methods 


Simulation of the JJA in the classical regime. At 
large u averaging out temporal fluctuations requires 
exponentially long times. Moreover, at u ^ 1 the 
resistance increases with u factorially leading to a 
strong heating in the computation of resistance un¬ 
less the measurement current is factorially small. Ob¬ 
servation of a small current against the background 
of a low frequency noise requires increasingly long 
times. Accordingly, for the realistic evolution times 
Tev ^ 10® the resistance can be computed only for 
u < 4.0. 

Simulation of the quantum problem. The time de¬ 
pendence of the entropy and the charge fluctuations 
for system of sizes L = 4,6,8 was computed using 
exact diagonalization in a symmetric subspace under 
charge conjugation. The tDMRG method was em¬ 
ployed for larger sizes which accuracy limits the range 
of times that we could study. In all simulations we 
impose the particle number conservation and open 
boundary conditions. 
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1 Free energy of the chain in 
thermodynamic equiiibrium 


Evaluation of the free energy with Hamiltonian § 
gives 

F=lT\n(^^'^-Tln[Io{Ej/T)] 


which allows us to express u in terms of T: 


u = 


T 


/ h{Ej/T) \ 
\ Io{Ej/T)) 


( 11 ) 


Here Ia{x) are modified Bessel functions. Using 0 
one can check the validity of <§ , which takes the form 
w = T'^du/dT . 


2 Noise distribution. 


At high temperatures the frequencies of individual 
phases are typically large, Wi S> 1 which allows to 
solve the equations ([^ in perturbation theory. Fur¬ 
thermore, because the effect of noise decreases expo¬ 
nentially with distance one can use the forward prop¬ 
agation approximation in which the evolution of the 
next phase is determined exclusively by the previous 
one 


dr"^ 


= -sin((()j - 4>i-i) 


( 12 ) 


Looking for the solution in the form (j)i = LOiT -f 
S(pi{T) and introducing the notation 


e 


i(pi _ 


/ 


duje 


for the noise at site i we get 

= 7-^Imz,_i(cc)e*(‘"-‘"')" (13) 

(wi -ojy 


We are interested in the propagation of the low fre¬ 
quency harmonics of the noise,cu ^ uii. Neglecting 
the high frequency components of the noise one ob¬ 
tains from (|13|) the noise recursion 


Zi{uj) 


2{uJi - w) 2 ' 


i-l(w) 


(14) 
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Iterations of (14) in the limit uj ^ uji lead to the equa¬ 


tion (5) of the main text. The exponential decrease of 
the noise away from its source in the classical regime 
is in agreement with the numerical computation in 
the quantum regime (see section]^. 

The recursion (5) of the main text implies that the 
noise generated by a single triad at distance, r 1 is 
given by a product of a large number of factors. Being 
a sum of a large number of random independent terms 
the logarithm of the noise has Gaussian distribution. 
Thus, at the distance >> r ;g> 1, from the closest 
triad the noise distribution has a log normal form 


PriO = 


A/2^rC(((r) 


exp 


(C-C*(r)f 

2 C*(r) 


(15) 



Fig. 8: Logarithm of the dimensionless amplitude of 
the charge transfer, /(r) determined in the 
forward propagation approximation. 


where C = Inx, C*(’') = (7 + Inu) r, and 7 = 0.577 is 
Euler constant. The exponential decay of the effect 
of a single triad ( 5|15 ) implies the charge localization. 
The DC charge transport in a macroscopically large 
array requires interaction of different triads through 
the quiet regions. A quiet interval (j,j -I- r) be¬ 
tween islands can be characterized by the resistance 
Rj,j+r ~ Convolution of distribution (15) 

with the Poisson distribution for the sizes of the quiet 
regions, r, yields the log-normal distribution for the 
resistances Rj,j+r of quiet regions. 


3 Charge propagation in quantum 
insuiating phase. 

In the insulating phase the charge transfer by dis¬ 
tance r appears in the order of the perturba¬ 
tion theory in the parameter Ej / ^/TEc- For small 
Ej j\/TEf. <C 1 the main contribution to the ampli¬ 
tude, 4'(r) of this processes comes from exactly r — 1 
virtual hops between neighboring sites in the same di¬ 
rection (forward propagation approximation). These 
hopes can occur in any order. Each sequence of the 
hops corresponds to a particular permutation, P, of 
1... r — 1. After n < r — 1 hops the charging energy 
changes by AE„ = E^ J2k<ni^P{k)+i-Qp{k)), so that 

(ln|vI/^^H^B(r)|) = 2E,\j:k{qPik)+i- mk))\) 


where ■ ■)g means average over charge configu¬ 
rations normally distributed with variance T/Eq- 
Rescaling the charge and factoring out the dimen- 
sionfull parameters we find that 

(ln|^FWD(r)|) =r-ln -H/(r) (16) 


fir) 



1 

2 \ J2ki<ip{k)+i - qp{k)) 


(17) 


where average is performed over charge configura¬ 
tions with unit variance. The numerical evaluation 
of /(r) shows that /(r) « ijr with 77 « —0.3(see Fig. 

©• 

We estimate the transition temperature by de¬ 
manding (In |4('i;'WD(^)|) ~ const which gives 

rpFWD ^ e^'^Ej/Ec- The forward propagation ap¬ 
proximation neglects the effects of the backscattering 
and level repulsion. All these apparently small ef¬ 
fects favor the insulating behavior, thus, the actual 
transition is below Tc ^ . We conclude that 

the numerical prefactor, rj, in the equation for the 
critical temperature, Tc = fjEj/Ec is close to one. 


7 ^ 0.5. 

The estimate of the charge transfer (16) neglects 
the charge discreteness which is expected to become 
irrelevant at T/Ec —>- 00 . To estimate the effect of 
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the discreteness at finite T/Ec we note that it is dom¬ 
inated by the appearance of exactly zero l^En for 
some n (resonances). Such degeneracies are lifted in 
the next order of the perturbation theory: the energy 
difference between the states with charge < 7 ^ -I- 1 and 
charge qi at island i becomes 

E{qi+ ^) — E{qi) = Ec{qi +-j) + 

Si = Si^i+i + 




g*+i - qi 


r)-| 2 


4 ( 9 *+! 9* 2 ) (9i-|-l 9i 2 ) 


Close to the transition line we estimate Si ~ 
{Ec/EjY, so the degeneracy is lifted by SE ~ 
E'q/Ej <C Ej. These resonances occur with the 
probability Ec/Ej, at an average distance r ~ 
EjjEc from each other. Close to the transition the 
hopping amplitude between these sites is given by 

«'(r) ~ iijexp (- 



Fig. 9: Time evolution in the insulator, Ej/Eq = 
0.05. The main panel shows slow evolution 
of the entropy which is roughly linear in Int 
in a wide range of times. The insert shows 
slow relaxation of the charge which becomes 
even slower for longer samples. 


4 Charge and entropy relaxation in 
quantum regime 


The energy difference, 
provided that 

T-Tr 


AE, exceeds this amplitude 


Tc ^ Ej Ec 


(18) 


where 7 ^ 1 is numerical coefficient. The many body 


localization takes place only when the condition (181 


is satisfied. We conclude that the shift of the transi¬ 


tion line upwards is small in Ec/Ej (18). 

According to 0 the charge propagation is con¬ 
trolled by the average logarithm of the charge dif¬ 
ference. This suggests that the simulations in which 
for the charge is evenly distributed in the interval 
{—Q, Q) should resemble the simulations for the nor¬ 
mal distribution with variance T/Eq when the two 
cases correspond to the same(ln( 7 i) or (ln(( 7 i — qi+i)) ■ 
Application of such criteria leads to the conclusion 
that restriction of the charges—Q < qi < Q) at 
T —)■ oot is equivalent to the unrestricted charges 0 
with effective temperature T = j'EcQ^ with 7 ' « 1 


Here we give the details of the numerical results for 
charge and entropy evolution in the insulating and 
metallic regimes that are both qualitatively different 
from the non-ergodic bad metal phase. 

We begin with the insulating phase. The behavior 
of the entropy in a wide time range {t < 10 ^°) is 
shown in Fig. one can distinguish three distinct 
regimes. 

At shortest times {Ect < 20) the entropy growth 
fast. The fact that this growth is identical for the sys¬ 
tems of different sizes suggests that it is due to the 
particles spreading across the boundary which hap¬ 
pens at time scales t ~ 1/Ej. In the intermediate 
time region (20 < Ect < 10 ^) the entropy is almost 
size independent for all but smallest sizes and small 
{S ^ 1). Similar behavior is observed for disordered 
systems; it is due to a small entanglement of typical 
excitations a bit further from the boundary or to a 
larger entanglement of rare excitations with excep¬ 
tionally close energies. Because the amplitude of the 
entanglement decreases exponentially fast with dis- 
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Fig. 10; Charge relaxation and entanglement en¬ 
tropy relaxation in a good metal realized at 
EjjEc = 1.75 


tance the contribution to the entropy coming from 
the states far from the boundary is small and thus 
the entropy demonstrates weak size dependence in 
this regime. 

The striking feature of the insulating state in this 
model is a very slow (logarithmic) growth of entropy 
at long times {Ect > 10^) that remains much below 
In 5 per site. We attribute this dynamics to degener¬ 
ate charge configurations: even when located far from 
each other such configurations can hybridize due to 
a small but non zero Ej/Eq- For example, config¬ 
urations connected by inversion symmetry possess a 
matrix element which is exponentially small in sys¬ 
tem size. The time needed to resolve this couple is 
thus exponentially large. 

In contrast to the entropy, the charge fluctua¬ 
tions display simple monotonic behavior at all times. 
The characteristic charge relaxation time increases 
extremely rapidly with the system size. We con¬ 
clude that this phase is a genuine insulator which 
phase space is separated into thermodynamically 
large number of independent compartments. 

In a good metal there is only one regime for charge 
and entropy time evolution. The charge relaxes 
quickly and this relaxation does not show any sign 
of getting slower at larger system sizes, see Fig. 
Accordingly, the entropy increases rapidly and sat¬ 
urates at the values that approach Lin5 for large 


Fig. 11 ; The insert shows the fidelity of the ground- 
state as a function of Ej/Ec for different 
sizes L. The minimum of the fidelity as 
a function of [ln(L)]“^ appear in the main 
panel. The red line fits the data to a line. 
The location of the critical point can be ex¬ 
tracted from the value of the fitting at origin 
Ej/Ec = 0.63 ± 0.04. 


system sizes. 


5 Berezinski-Kosterlitz-Thouless 
critical point at zero temperature 

In this section we give the details of the numerical 
methods that allowed us to determine the value of 
the ratio of Josephson to charging energies, t], for 
the Berezinsky-Kosterlitz-Thouless phase transition 
of the model Q at zero temperature. Finite order 
derivatives of the energy do not display any discon¬ 
tinuity at this transition. Thus, we have to use more 
sophisticated numerical methods than the ones used 
in the case of second order phase transitions. 

The idea of our approach is to compute the fidelity 
of the ground-state defined by the equation 


E{Ej/E,) =< ^gs{Ej/E, + SMgsiEj/E,) > (19) 

for small <5. The result obtained by DMRG method 
is shown in the insert of Fig. 11 as a function of 
Ej/Ec for different sizes. One expects that fidelity is 
minimal at the critical point |26j . In the main panel 
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of Fig. 11 we show positions of this minima as well as 
its extrapolation to infinite size using Ej/Ec{L) = 
a\ln{L)]~'^ +7? [27]. This procedure yields r] = 0.63 ± 
0.04 cited in the main text. 




